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Mapping complex problems to simpler effective models is a key tool in theoretical physics. One 
important example in the realm of strongly correlated fermionic systems is the mapping of the 
Hubbard model to a t-J model which is appropriate for the treatment of doped Mott insulators. 
Charge fluctuations across the charge gap are eliminated. So far the derivation of the t-J model 
is only known at half-filling or in its immediate vicinity. Here we present the necessary conceptual 
advancement to treat finite doping. The results for the ensuing coupling constants are presented. 
Technically, the extended derivation relies on self-similar continuous unitary transformations (sCUT) 
and normal-ordering relative to a doped reference ensemble. The range of applicability of the 
derivation of t-J model is determined as function of the doping S and the ratio bandwidth W over 
interaction U. 

PACS numbers: 71.10.Fd, 75.10.Jm, 71.27.+a, and 71.30.+h 



I. INTRODUCTION 

The Hubbard modeP^ is one of the most common 
models for the description of strongly correlated electron 
systems on lattices. Because it contains the motion of the 
electrons as well as the interaction between two electrons 
at the same site it is capable to describe charge degrees 
of freedom as well as magnetic degrees of freedom. Due 
to the rich physical behavior of the Hubbard model an 
analytic solution is not possible except in one dimension^. 

One common route to simplify the model for large re- 
pulsion U is to derive an effective model which does no 
longer contain charge fluctuations across the charge gap. 
Processes which change the number of doubly occupied 
sites (double occupancies, DOs) are eliminated. For large 
enough repulsion U and at half-filling the electrons are 
fixed on their lattice sites. In this Mott-insulating phase 
the model can be mapped onto a Heisenberg model de- 
scribing only the energetically low-lying spin degrees of 
freedom. In the immediate vicinity of half-filling the mo- 
tion and the interaction of doped holes is described by the 
extension of the Heisenberg model to the t-J modePH^l. 
The metallic behavior for small values of the repulsion is 
beyond the applicability of this mapping^. 

In the present work the mapping of the Hubbard model 
to the t-J model is extended to finite macroscopic dop- 
ing concentration 6. The influence of the doping 6 on 
the resulting parameters of the t-J model is studied for 
sufficiently large repulsion. Our approach provides a sys- 
tematic and controlled derivation of the effective cou- 
pling constants as function of the doping concentration. 
Thereby, an important gap between the applicability of 
the derivation of the t-J model and its actual applications 
is closed. 

First, we consider the half-filled case. The elimination 
of the charge fluctuations across the charge gap is per- 
formed by a self-similar continuous unitary transforma- 
tion with various types of generators. Besides the mag- 
netic exchange couplings, the resulting t-J model con- 



tains the motion and the interaction of holes and doubly 
occupied sites. Since the mapping starts from a refer- 
ence ensemble comprising the two spin states with equal 
weight and without any correlations between the spins 
on neighboring sites the spin state in the effective model 
remains unspecified. 

The mapping relies on the elimination of processes 
changing the number of holes and doubly occupied sites. 
Note that an empty site represents a double occupancy 
of two holes. In the half-filled case the density of states 
exhibits two distinct bands for large U (see Fig. [IJ . The 
bands display eq ual w eight 1/2 and they are well sepa- 
rated for large £fl2UlHSI Thus the states without holes 
or doubly occupied sites are energetically well separated 
from the ones with one or more holes or doubly occu- 
pied sites. If U is decreased the bands approach each 
other. As soon as they touch the insulating phase is no 
longer the appropriate phase and metallic behavior oc- 
curs resulting in the breakdown of the mapping to the 
t-J modelP 
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FIG. 1: (Color online) Density of states for the half-filled Hub- 
bard model with large repulsion L^"^. The density of states 
exhibits two distinct, equally weighted bands, the low er Hu b- 
bard band (LHB) and the upper Hubbard band (UHB jUnES 

The generic density of states obtained in the doped 
case is depicted in Fig. [2} The effect of the doping on the 
density of states consists in shifting the Fermi energy into 
the lower Hubbard band for hole doping and redistribut- 
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FIG. 2: (Color online) Density of states for the case of hole 
doping with a doping concentration S. The weight of the hole 
state is given by 8. The two half-filled states both carry the 
weight . Note that the weight of the lower band is larger 
than the one of the upper Hubbard bancP. 



ing the weight of the bands. Electron doping is com- 
pletely analogous in shifting the Fermi energy into the 
upper Hubbard band. Since we focus here on particle- 
hole symmetric, bipartite lattices we will consider hole 
doping without loss of generality. 

It is obvious from the comparison of Fig. [T] with Fig. 
[2] that the energetic separation of the Hubbard bands is 
more subtle in the doped case than in the half-filled case. 
The bands are shifted depending on 5 and spectral weight 
is transferred as well. The shift of spectral weight is a 
smoking gun evidence for strongly correlated fermionic 
systems. 

Simple counting arguments in the limit U — > 00 tell 
us the distribution of weight. Adding an f electron to 
a site succeeds with probability p = 5 + (1 — 5)/2. The 
first term results from the fraction of empty sites; the 
corresponding weight is found at low energy because no 
doubly occupied site has to be created. The second term 
results from the fraction of sites occupied by J, electrons; 
the corresponding weight is found at about lu s» U be- 
cause a doubly occupied site is created. Together with 
the sum rule the weights shown in Fig. [2] result. 

We continue to use the number of doubly occupied sites 
(DOs) as criterion to distinguish different sectors of the 
Hilbert space. In order to have a quantitative measure for 
the energy separation of the sectors with differing number 
of DOs the apparent charge gap A g is introduced which 
measures the energy separation of subspaces. It does not 
measure the energy gap between two pure states which 
is the reason why we call this separation of energy scales 
"apparent". While the apparent charge gap is not an 
energy gap in a rigorous sense it is experimentally signif- 
icant: It quantifies the energy needed in an Mott insula- 
tor to create a charge excitation irrespective of the spin 
state of the system. For instance, the system can be at 
a temperature which implies a disordered paramagnetic 
spin state while it preserves the insulating properties. 

In the half-filled case the apparent charge gap measures 
the minimal energy of a doubly occupied site moving in 
an arbitrary spin background. The apparent charge gap 
is reduced for increasing values of the band width W. If 



the gap vanishes the system is no longer stable against 
charge fluctuations and the mapping fails^. 

Obviously, the physical properties of strongly interact- 
ing fermionic systems depend considerably on the doping 
levePS. This fact leads automatically to the question how 
the validity of the mapping from the Hubbard model to 
a generalized t-J model is influenced by doping. Thus 
the apparent charge gap has to be determined in depen- 
dence on the doping 5. Keeping track of the apparent 
charge gap we determine the parameter range in which 
the mapping is still justified. This range of applicability 
may not be misinterpreted as phase diagram although it 
has some similarities^. For instance, for large repulsion 
U a doped t-J model is still perfectly well defined while it 
displays metallic behavior. Of course, it is expected that 
the applicability of a t-J model decreases upon increasing 
doping^!. 

In view of the above, it is one of our central objectives 
to derive a diagram showing the range of applicability in 
dependence on the doping which has, to our knowledge, 
not been done before. Our findings provide access to 
the limitations of the use of t-J models in the context of 
planar cuprates to the extent that they can be described 
by a single-band Hubbard model. 

The approach used is based on two conceptual ingredi- 
ents. The first is a systematically controlled change of ba- 
sis by means of continuous unitary transformations^ 9 -^. 
The second is the choice of a doped reference ensemble 
without spin or charge order. Within the range of appli- 
cability the effective t-J model is derived. The doping 
dependence of the effective coupling constants is studied. 
The results are given in dependence on the ratio W/U 
and on the dopant concentration 8. The method im- 
plemented here uses a self-similar truncation scheme to 
reduce the amount of proliferating terms in the running 
Hamiltonian. The truncation is performed according to 
the range of the processes, that means rather local pro- 
cesses are kept while ones of longer range are neglected. 
Hence the local processes acquire a non-perturbative de- 
pendence on the bare, initial coupling constants of the 
system. 

Furthermore, recently introduced modified generators 
of the CUT are implemented to cope with the vast 
amount of terms arising during the transformation^. 
Their results are very clos e to t he previously used 
particle-conserving generato r^ * 18 * 19 !, but they signifi- 
cantly facilitate the calculation in terms of required mem- 
ory and CPU time. 

After this introduction, the model and the method (see 
Sects. |TTJ[TTT) are introduced. In Sect. |VI| results for the 
apparent charge gap are presented and Sect . | VII] provides 
exemplary results for the doping dependence of the cou- 
pling constants. Sect. VIII concludes the article. 
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II. HUBBARD MODEL 

We consider the fermionic Hubbard modeP^. It de- 
scribes electrons with spin <j on a lattice site i by their 
creation operator c\ a and their annihilation operator 

c ia . The Hamiltonian consists of two terms describ- 
ing the single-fermion kinetics (H t ) and their interaction 
(Hu) 



H = H t + Hu 
Ht = tJ2 (4sv + h-c) 



<»j> 



(la) 
(lb) 

(lc) 



The kinetic part consists of the hopping of an electron 
with spin a from site i to site j and vice versa. For 
this process to take place i and j have to be nearest 
neighbors as indicated by the bracket under the sum. 
The corresponding matrix element is denoted by t. The 
band width W of the model is given by W = 2zt with the 
coordination number z (number of nearest neighbors). 
In this work the lattice studied is the two dimensional 
square lattice with coordination number z = 4 so that 
W = 8t. 

The second part of the Hamiltonian determines the 
interaction of the electrons. This term constitutes a pure 
on-site interaction. In Hu the operator n,- )CT = c\ (T c i a 
represents the number operator for the electrons. This 
indicates that putting two electrons on the same site costs 
the additional energy U. 

In the Hubbard model there are four possible states per 
site. The site may be singly occupied by one electron with 
spin up or spin down | l), doubly occupied by two 
electrons with opposite spin | or completely empty 
|0). The last two configurations correspond to charge 
fluctuations and are referred to as double occupancies 
(DO) in this context. 

The interplay of motion and interaction of electrons in 
the Hubbard model provides a description of the metal- 
insulator transition 20 . Another important field of appli- 
cation of the single -ban d Hubbard model is the physics 
of high-Tc cupratefP^l. 

For a large Hubbard repulsion U in the half-filled case 
the density of states exhibits two separate bands, see 
Fig. [I] the so-called lower (LHB) and the upper Hub- 
bard band (UHB). For infinite U each site of the lattice 
is occupied by one electron which is energetically fixed to 
its site. If U is finite the electron can move and virtually 
hop to adjacent site. Thereby, DOs are created but the 
the physics remains rather local, that means, the charge 
correlation length stays small. 

Based on the locality of the important processes one 
may map the Hubbard model onto an effective t-J model. 
The generalized t-J model conserves the number of DOs. 
It comprises a part which describes the magnetic degrees 
of freedom, which is a generalized Heisenberg model, and 



a part which describes the motion and interaction of DOs 
reflecting the charge degrees of freedom. In order to ob- 
tain a model conserving the number of DOs, processes 
which create or annihilate DOs have to be eliminated. 
One systematic way to achieve this objective is the ap- 
plication of continuous unitary transformations to the 
Hamiltonian. 

For smaller values of U the local picture used here is no 
longer appropriate and the derivation of the t-J model is 
not justified. 



III. CONTINUOUS UNITARY 
TRANSFORMATIONS 

A. General Framework 

The effective t-J model is derived from the Hubbard 
model by eliminating processes which change the number 
of DOs. This elimination is pe rform ed using continuous 
unitary transformations (CUTpl2USHl2l xhe elimination 
is based on a systematic change of the basis 



H{£) = U (£)H U\l) 



(2) 



with a unitary operator U and a continuous auxiliary 
variable £ referred to as the flow parameter. The trans- 
formation is determined by the flow equation 



H(£) = [ v (e),H(e)] 



(3) 



where rj{£) denotes an antihermitian infinitesimal gener- 
ator. At I = the transformation starts with the ini- 
tial Hamiltonian H. The unitary transformation can be 
stopped at any arbitrary value of the flow parameter £. 
Usually, the effective Hamiltonian is reached for £ = oo. 
Due to the continuity of the transformation it is read- 
justed to the flowing Hamiltonian for every value of £. 

The transformation stops automatically when the com- 
mutator [H(£), r/(£)) vanishes which is generically the case 
for £ — )■ oo, i.e., for convergence for £ — > oo. The struc- 
ture of the effective Hamiltonian is determined by the 
choice of the generator r](£). We first choose the gen- 
erator which leads to an effective model conserving the 
number of DOs. To this end, we introduce the operator 

d ■= E[™a™a + 0- - ™a) (! - "a)] ( 4 ) 

counting the number of DOs. 

By the use of D the repulsive part of the Hamiltonian 
can be written as 



(5) 



with N denoting the number of sites. The kinetic part 
is split into three parts according to their effect on the 
number of DOs 



Ht — Tq + 7+2 + T_2 



(6) 
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where Tj creates i DOs. The terms are given by 



<l,J>,cr 



(7a) 



<i,j>,er 



'J,°4,<A<f (1 ~ (7b) 



= *-5 



1 - ^t.o-) c\ s C j ^flj,a + 



<i,j>,o 



(7c) 



with a := —(7. 

The terms contained in To have no effect on the number 
of DOs while the terms in T + 2(T_2) increase (decrease) 
the number of DOs by two. These terms are the ones 
that we intend to eliminate by the transformation. In 
the initial Hamiltonian the prefactors t + 2,t— 2 and to are 
equal, but they evolve differently under the CUT. 

The first generator we use, the so-called quasiparticle 
conserving generator r] pci can be expressed by the com- 
mutator 



Vpc (£)= D,H(£) 



(8) 



This generator corresponds to the generator defined in 
Refs. l9H7fHT)l except for a global factor of two, which 
just implies a multiplicative renormalization of the flow 
parameter. The terms comprised by this generator are 
sketched in Fig. [3] The terms of the Hamiltonian are 
classified according to their number of quasiparticle cre- 
ation and annihilation operators. The (J, I) block consists 
of terms with j creation and I annihilation operators. 
Note that such a term requires at least I excitations to 
be present in order to become active. But it is also active 
if more than I excitations are present in the system. In 
this respect, the scheme in Fig. [3] may not be mistaken to 
be a matrix. For a comprehensive presentation we refer 
the reader to Ref. [T71 



The quasiparticle conserving generator rj pc comprises 
all terms of the off-diagonal blocks. Due to the structure 
of the generator the block-band str ucture o f the Hamil- 
tonian is preserved during the flo W 18 * 19 ^. During the 
whole flow there will only be terms created which change 
the number of DOs by 0, +2 or —2. 

With the definition Q of the generator the flow equa- 
tion ([3| can be calculated. Comparing the contributions 
on both sides of Eq.[3]a set of differential equations for the 
prefactors of the monomials in the creation and annihi- 
lation operators is obtained. These differential equations 
are first order in £ and they are bilinear in the prefactors 
entering on the right hand side. 

The equations do not form a closed set because in in- 
finite systems new terms continue to arise on each appli- 
cation of the commutator in Q. For i = these terms 
carry the prefactor zero because they are not part of the 
initial Hamiltonian. If we kept all these new terms in the 
remaining calculations we would obtain exact results for 
the effective model. But the number of arising terms is 
rising exponentially so that we have to limit them in num- 
ber. For this purpose a truncation scheme is introduced 
which specifies the relevance of all term. Less important 
terms arc neglected, leading to a closed set of differential 
equations which can be solved numerically. 

In many previous applications a small parameter is 
used to classify the arising terms so that a perturbative 
treatment results. In contrast, we are adopting here a 
truncation scheme which classifies the terms according 
to their structure. Such a CUT scheme is usually called 
self-similar. It resembles more conventional renormaliza- 
tions. Effects of infinite order are present in the prefac- 
tors of the kept terms. 

The truncation scheme used in this work keeps or ne- 
glects terms according to their locality, i.e., according to 
the range of the represented physical process. This ap- 
proach is well justified if the model under study is gov- 
erned by a small correlation length. This is exactly the 
case for a Hubbard model at large U where the propa- 
gation of charge degrees of freedom is suppressed by the 
high energetic cost of creating a DO. 




FIG. 3: (Color online) Schematic diagram of the Hamiltonian 
with the terms contained in the n vc generator highlighted in 
red (dark grey). These terms will be eliminated. The terms 
in the uncoloured squares are zero initially and stay zero in 
the flow induced by n pc . 



B. Reference Ensemble and Normal Order 

Before we discuss how we measure the degree of locality 
we have to find a unique representation for the operators 
to be sure to treat similar terms in the same way. To 
this end, the monomials are expressed as normal-ordered 
products of local operators. The normal-ordering we are 
using is not the standard one known for the fermionic or 
bosonic algebra because the creation or annihilation of a 
DO can not be represented by interaction free fermions 
or bosons. Instead, we use a reference ensemble. Non- 
trivial operators are only those which create or annihilate 
fluctuations away from the reference ensemble. For a 
given doping concentration 5 the reference ensemble is 
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bosonic 


fermionic 


1 

cr z = nf — hi 
c 1 c l 

hs = n-f + hi — 1 + (51 
Ds = 2h-fhi — fig 


(1 - hi) c t 
(1 - hi) &i 

hih 

hlCi 

hic\ 

%cj 
(1 ~ ni) cj 
(1 - fi t ) c{ 



TABLE I: Basis of normal-ordered local operators 



denned by the statistical operator 
1-5 r 



ps = 



tMtii+uMH +mi(o\i> (9) 



where the product extends over all lattice sites i. In the 
half-filled case (5 — 0) the reference ensemble is param- 
agnetic and the magnetic degrees of freedom are totally 
disordered. Each site is equally probably occupied by 
an j" or by a \. electron; no direction is singled out, no 
correlation between neighboring sites exists. Charge fluc- 
tuations from this reference ensemble are the empty |0) 
and the doubly occupied site | J/f). Magnetic fluctua- 
tions are induced by the application by spin operators, 
see below. 

Considering doping we focus on hole doping only be- 
cause the model at hand is particle-hole symmetric so 
that electron doping leads exactly to the same results. 
Hence we include the empty state |0) in the reference 
ensemble ^ besides the half-filled states with a proba- 
bility given by the doping level 5. The remaining weight 
is again equally distributed over the two spin states. Note 
that this extension to the doped case does not introduce 
any bias. There is no correlation between sites nor on 
each site. Hence the reference ensemble ^ is the mix- 
ture with the maximum entropy at given level of doping. 

Based on the reference ensemble we define a term as 
normal-ordered if the expectation value of each of its fac- 
tors of local operators vanishes with respect to this en- 
semble. Thus a normal-ordered local operator fulfills 



= 0. 
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(10) 



Based on this condition we define a basis of local normal- 
ordered operators (see Tab. [I]). Of course, the identity 
does not fulfull the condition ( |To| . But the identity is 
the trivial action of an operator and obviously does not 
create or annihilate any fluctuation away from the refer- 
ence ensemble. Without the identity the list of operators 
would not be complete. Any monomial occuring during 
the flow is expressed in the operator basis in Tab. [I] 

Among the normal-ordered operator basis the operator 
hs = h\ + hi — 1 + 51 occurs. This operator counts the 
number of electrons on one site relative to the mean value 



of the filling 1 — 5. In the half-filled case (5 = 0) the mean 
value of the filling is 1. Thus ho applied to an empty site 
yields —1. Applied to a doubly occupied site yields +1 
and a singly occupied site leads to 0. In the doped case 
the counting operator Ds can be determined from the 
operator for the half-filled case through 



D s = D a + Sj2^ 



(11) 



A unique representation for a possible operator occur- 
ring in the Hamiltonian or in the generator is given by 
the appropriate linear combination of monomials of the 
basis operators. The monomial is the product of local op- 
erators acting on different sites. Hence the expectation 
value of each monomial also vanishes. 



C. Implementation: Truncation Schemes 

The cluster of sites of a monomial is the set of sites 
on which the monomial has a non-trivial action. Here 
'non-trivial' simply means no to be the identity. Based 
on the clusters a truncation scheme is defined by mea- 
suring its locality by the extension of its cluster. The 
extension is defined as the maximum taxi cab distance 
between the outermost cluster sites. Thus the extension 
in x and in y direction have to be summed up. An ex- 
emplary term with an extension of 3 is shown in Fig. 
[4] With the normal-ordered operators given in Tab. [TJ 
a monomial with the cluster shown in Fig. [4] can be ex- 
pressed as product of 3 local operators. These operators 
act on the lattice sites (0,0), (2,0) and (2,1). The exten- 
sion of this cluster in x-direction is 2 and its extension in 
y-direction is 1. 



FIG. 4: (Color online) Cluster of a term with an extension 
of 3. The term consists of the operators A,B and C The 
taxi cab distance between the outermost operators A and C 
determines its extension. 

To limit the number of generated terms in the course 
of the flow we define a maximum extension. For each 
normal-ordered term generated by the commutator the 
extension is determined. A term with an extension higher 
than the defined maximum extension is neglected. A 
CUT truncated to a maximum extension of two only con- 
siders terms whose clusters have an extension two or less, 
(see Fig. [sj. 

In this way more and more extended truncation 
schemes are used until the numerical results do not 
change noticeably anymore. Then the calculation is slop- 
pily said to be 'converged'. To illustrate how the cou- 
plings change under the influence of different truncation 
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FIG. 5: (Color online) The maximum clusters occurring in a 
calculation with a maximum extension 2. In the following we 
will call the calculation based on this truncation the plaquette 
calculation. 

schemes we consider the nearest-neighbor magnetic ex- 
change constant J± 

^Heisenberg = J\ ^ ' S%Sj. (12) 

d,i) 

In leading perturbation order one obtains j[ 2 ^ = ^JHEI. 
The results for this coupling constant obtained by CUT 
are shown in Fig. [6j 



sufficient to describe 6th order processes of J\ . Since the 
double plaquette calculation results in a large number 
of terms an additional truncation scheme is introduced, 
the l upto4' truncation. In this calculation a subset of 
processes with extension 3 are considered which consist 
at most of 4 non-trivial local operators on different sites. 

From the results for J\ we deduce that J\ already con- 
verged for a maximum extension 2. Thus the terms con- 
tained in a plaquette calculation are sufficient to describe 
the nearest-neighbor exchange appropriately. Consider- 
ing even more extended terms does not change the result 
significantly. 

We emphasize that the calculation with a truncation 
according to the extension is not equivalent to a finite- 
size cluster calculation. The latter acts on a finite cluster 
with finite Hilbert space only. The former only restricts 
the maximum range of physical processes but remains a 
calculation on the thermodynamic, infinitely large sys- 
tem. The latter computes quantities on finite clusters 
and a second approximation, for instance finite-size scal- 
ing, is needed to extend these finite-cluster results to the 
infinite system. 



D. Implementation: Flow Equations 




0.5 1 1.5 

w/u 



FIG. 6: (Color online) Effective nearest-neighbor Heisenberg 
exchange Ji obtained in various truncation schemes relative 
to the leading perturbative result in t/U. 

The results are shown for various truncation schemes, 
where 'min' denotes the minimal model in which only 
the Heisenberg exchange term is kept in addition to the 
terms present in the initial Hamiltonian. The NN trun- 
cation represents a nearest-neighbor calculation defined 
by the maximum extension 1. This calculation repro- 
duces the second order perturbative result for J\. The 
plaquette calculation contains all terms which fit on the 
clusters shown in Fig. [5] This truncation corresponds to 
a maximum extension 2. It reproduces 3\ up to fourth 
order in t/U. 

A maximum extension of 3 corresponds to the so-called 
'double plaquette' calculation. This truncation scheme is 



The choice of a truncation scheme implies that the 
set of differential equations describing the flow is finite. 
Hence, two tasks have to be accomplished, both of which 
are implemented on computers. First, the flow equations 
have to be set up. Even though the large number of run- 
ning coupling constants makes the use of computer aid 
indispensable, this step is an essentially analytic calcula- 
tion. Second, the flow equations are integrated numeri- 
cally which results in the effective coupling constants de- 
termining the effective model. In the effective model the 
most important subspaces, subspaces of different number 
of DOs, are decoupled from the rest of the Hilbert space. 
Thus important observables can be calculated with less 
effort. 

The derivation of the flow equation is realized by a 
program implemented in CH — h This program performs 
the calculation of the commutators and collects all con- 
tributions to the same term. Due to the vast amount of 
terms it is advantageous to use symmetries to increase the 
efficiency. If a particular term can be generated from an- 
other term by applying symmetry transformations both 
terms have the same pref actor. The model under study 
displays the SU(2) spin rotation symmetry and the point 
group symmetry of the square lattice. This group con- 
tains rotation symmetries about 7r/2, 7r and 3/2-7T, re- 
flection symmetries about x, y and the diagonal. In the 
half-filled case the particle-hole symmetry may be used 
additionally. Of the spin rotation symmetry we only ex- 
ploited the spin flip symmetry, i.e., the U(l) symmetry 
of rotations around S z . In addition, we used that the 
Hamiltonian is hermitian conjugate so that adjoint terms 
also must have the same prefactor. 
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By applying the above symmetries up to 64 terms are 
created out of a single term. Since they all carry the 
same prefactor it is sufficient to treat one representative 
instead of all 64 terms separately. By this technique, the 
number of terms is reduced from more than 1.6 million 
to 26251 in the double plaquette calculation. Yet the 
double plaquette calculation remains costly. It requires 
14.7 weeks of CPU time and more than 20 GB RAM 
memory. 

Compared to the derivation of the flow equations the 
solution of them is straightforward. We start at I = 
with the initial Hamiltonian and integrate the differential 
equations. At t = oo the effective model is reached. Since 
the integration is performed numerically, this limit can 
not be reached and we stop before at large enough values 
of I. 

In order to have a measure to which extent the CUT 
is accomplished we introduce the residual off-diagonality 
(ROD) 13. The name is motivated by the idea that the 
CUT eliminates the off-diagonal terms. As we will see 
in the next subsection the precise choice which terms are 
eliminated and which are not depends on the choice of 
the generator r\. Hence in practice the ROD is a measure 
of the norm of the generator. The ROD is calculated by 
squaring the (real) prefactors of the terms of the genera- 
tor, summing them and finally taking the square root of 
this sum. 

The ROD measures to what extent the terms in the 
generator are eliminated at the current value of the flow 
parameter I. When the ROD vanishes, the generator 
vanishes and consequently the transformation is finished. 
When the ROD is decreased to some small value, for 
instance 10~ 15 , the calculation can be stopped at t < oo. 
The contributions of the remaining off-diagonal terms are 
negligible so that we consider the model obtained to be 
the wanted effective model. 
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E. Various Choices of the Generator 



Because the number of generated terms during the flow 
leads to computational costly calculations, we consider 
various choices of generators for simplification^. The 
basic idea of the modified generators is that the most 
relevant physics requires only a very small number of 
DOs. Hence it may be sufficient to separate subspaces 
with zero or one DO from the remaining Hilbert space in- 
stead of applying rj pc which eliminates all terms changing 
the number of DOs. 

An obvious example is the derivation of the Heisen- 
berg model describing the magnetic degrees of freedom 
without any charges. Here the separation of the sub- 
space without any DO from the remaining Hilbert space 
is completely sufficient. Thus we consider the generator 



N 

%s W = 2£(i?^)-#°W) (13) 

i>0 



where N denotes the number of quasiparticles. The op- 
erator Hj (£) represents all terms which contain j annihi- 
lation operators of DOs and i creation operators of DOs. 
This generator contains all terms which couple to the 
subspace without DOs. Note that this subspace is a high- 
dimensional subspace and not a single ground state for 
the model under study in contrast to the situation con- 
sidered by Fischer et alP^. But the other conceptual 
points, e.g., concerning the formulation in second quan- 
tization and the differences to a matrix formulation^ are 
the same. The Hamiltonian and its evolution under the 



CUT induced by the gs-generator ( 13 ) is graphically rep- 
resented in Fig. [7j 
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FIG. 7: (Color online) Effect of the gs-generator. The terms 
of the Hamiltonian are labelled according to the number of 
creation and annihilation operators they contain. Thus a term 
in the {i, j} block creates i DOs after annihilating j DOs. 

If in addition we aim at an explicit description of the 
motion of a single DO the generator rj gs ^i p has to be used 



(14) 



jv 



N 



i>0 



i>0 



In this generator the idea of decoupling some subspaces 
from the remainder of the Hilbert space is extended to 
the subspace with one DO. Thus also terms coupling to 
this subspace are included as depicted in Fig.[HJ Whereas 
the subspaces with zero and with one DO are decoupled 
at the end of the transformation, the other subspaces are 
still coupled. To compute eigenvalues in the subspaces of 
zero or one DO only these subspaces need to be taken into 
account. In contrast, eigenstates involving two or more 
DOs still require the diagonalization of the full Hilbert 
space. 

The CPU time needed for a double plaquette calcula- 
tion using various generators are given in Tab. |Hj In the 
case of a CUT based on the gs,lp-generator the use of 
symmetries is even more efficient. Using all of them ex- 
cept for the particle-hole symmetry reduces the number 
of terms from 5 million to 55049. The reader may be sur- 
prised that these numbers are larger than those for the 
particle-conserving ?y pc although terms linking subspaces 
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FIG. 8: (Color online) Effect of the gs,lp-generator. The 
terms of the Hamiltonian are labelled according to the number 
of creation and annihilation operators they contain. Thus a 
term in the block creates i DOs after annihilating j 

DOs. 

pc-generator I gs, lp-generator I gs-generator 



102 days 



51 days less than 10 days 



TABLE II: Comparison of CPU time needed for a double 
plaquette calculation using all symmetries in the half-filled 
case with different generators. 



with higher number of DOs are not decoupled. The ex- 
planation is that ri-p C has the additional feature that it 
preserves the block-banded structure of the Hamiltonian 
while the other generators do nolPfflH. Yet the modified 
generators induce a simpler and faster CUT as shown by 
the numbers in Tab. [Til 



IV. MINIMAL AND NEAREST NEIGHBOR 
MODEL 

In this section we present analytic solutions of the flow 
equations which are possible for the two simplcmost trun- 
cation schemes. Due to the simplicity of the truncation 
schemes no difference between the different choices of the 
generator are found. For concreteness, we consider the 
particle-conserving r/ pc here. 



A. Minimal Model 

The calculation of the minimal model starts by study- 
ing all processes on between adjacent sites. This nearest 
neighbor (NN) calculation is equivalent to a maximal ex- 
tension of e = 1. To arrive at the minimal model all 
terms not present in the initial Hamiltonian except the 
NN Heisenberg exchange are neglected. 

The initial generator takes the form 



r,{£)= D,H(£) 



= 27! 



+2 



2T_ 



with the flow parameter £. Inserting this definition in the 
flow equation ^ we calculate 



d£ 



H(£) = r,(£),Hu(£) + H t (£) 



(16) 



From the commutator new terms arise^. In the minimal 
model all terms except the NN Heisenberg exchange 



H] 



Hciscnb. 



&i&3 = 2 ( a 



+ °7 



(hi) 
4 -i 



(17a) 
(17b) 



are omitted. The coupling constant Ji(£) starts at 
Ji(0) = because it is not part of the initial Hamil- 
tonian H(0). For £ 7^ it evolves according to the flow 
equation. 

In this simple case the flow equation can be solved an- 
alytically for a general lattice with coordination number 
#3. For £ = oo the effective model is reached. The 
effective coupling constants take the form 



£n,eff — to 

t+2,cff = 

1 

2 
1 

z 



t/cff 
Jl.eS 



4U$ 



IQztl 



'-u . 



(18a) 
(18b) 

(18c) 
(18d) 



To obtain the equations for the square lattice z = 4 must 
be inserted. The variables t and Uq represent the initial, 
unrenormalized values of the hopping and the Hubbard 
repulsion. For simplicity we will omit the subscript and 
label the unrenormalized values by t and U henceforth. 
Since the terms T +2 and T_ 2 are hermitian conjugates 
and we assume their coefficients to be real t + 2 = t-2 
holds. From (18b) we see that the terms contained in 
T + 2 and T_2 arc eliminated as it should be because they 
change the number of DOs. The effective model is even- 
tually given by 



1 ; 

H c ft — U c ff-D 



To + Ji 



,eff 

<i,j> 



S{Sj 



(19) 



B. Nearest Neighbor Model 

In the nearest neighbor model all terms arising from a 
nearest neighbor calculation are included in the effective 
model as well as in the generators. This is the full calcula- 
tion with extension 1. It can still be solved analytically^. 
In this truncation scheme the Heisenberg exchange cou- 
pling is given by 



J 



l.off 



3 + z 



v /[/ 2 +4(3 + z)^-[/ j 



(20) 



Note the differences to the result of the minimal model 

(15) 

(jl8d|) . Of course, this differences arises only beyond lead- 



ing order. 

Besides the Heisenberg exchange the calculation con- 
tains the term Hy describing the interaction of two DOs 



H v {£) = V{£) J2 n i n r 



(21) 



<i,3> 
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The operator n counts the amount of electrons compared 
to the mean value of the filling. For the half- filled case 
the mean value is l.The third term created during the 
flow is 



(22) 



<hj> 



describing pair hopping processes of DOs. One of the 
processes contained in this term is the hopping of two 
electrons from site j to an empty site i. As the empty 
state as well as the doubly occupied state represents a 
DO this process does not change the number of DOs. In 
the effective model V and V„ take the values 



V cS 

Vp^eff 



3 + 42 
4 



y/U2 + 4(3 + z)t2 + 



3 + 42 



y/U2 + 4(3 + z)t2 + 



3 + 4z 
4 



42 



U (23a) 
U . (23b) 



In the case of doping another contribution to the How 
equation also arises. It reads 



(24) 



and determines the chemical potential /i. In the effective 
model this constant takes the value 



Sz 

2(3 + 2)' 



Mcff= n/0 , _ A E/\/l+4(3 + 2)— - o/o ~ - US. (25) 



i 2 _ 2 
IP ~ 2(3 + 2)' 



with the coordination number z. In leading order in jj 
this yields a chemical potential which depends linearly on 
the doping constant 8 and on the coordination number 
of the lattice 



( 2 ) = SZ 



L 

U 



(26) 



V. 



INFLUENCE OF THE CHOICE OF 
GENERATOR 



In this section the influence of the choice of the gen- 
erator is studied. First, we consider the ROD defined in 



Sect. Ill for the gs-generator (13 1, gs,lp-generator (14 1, 



and the pc-generator ([8|). The results for the ROD ob- 
tained in a double plaquette calculation at half-filling are 
shown in Fig. [9] as functions of the continuous flow pa- 
rameter I. 

For the gs-generator, the ROD converges for all values 
of the ratio W/U. In contrast, the ROD for the gs,lp- 
generator and even more pronounced the ROD for the pc- 
generator show non-monotonic behavior for larger values 
of W/U. 

Non-monotonic behavior of the ROD suggests that the 
intended transformation does not succeed. There is no 
strict statement that a successful CUT has to have a 
monotonic ROD. It is well possible that the ROD displays 
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FIG. 9: (Color online) Behavior of the residual off-diagonality 
(ROD) in a double plaquette calculation for various choices 
of the generator for various initial ratios W/U. 



local maxima which indicate that some energy eigen- 
states are re-ordered, see for instance Ref. If the 
CUT is performed without approximation any unitary 
transformation is as good as any other. But since we 
have to truncate many terms the upturn of the ROD in- 
dicates a potential loss of accuracy. If the total norm of 
the off-diagonal terms is large there is still a significant 
transformation to be done. In the course of this transfor- 
mation the truncation of terms may introduce significant 
errors. In return, a quickly decreasing ROD indicates 
that all coefficients to be eliminated decay fast and sig- 
nificant truncation errors are less likely. But we like to 
stress that the behavior of the ROD is only an indica- 
tor for possible truncation errors which eventually may 
imply that the intended mapping breaks down. 

The faster convergence of the gs-generator is straight- 
forward to understand because the gs-generator com- 
prises only terms which create DOs from the reference 



ensemble or which annihilate them, see Eq. ( 13 ) and Fig. 
[7] As long as there is a finite charge gap A g these pro- 
cesses are exponentially suppressed: cx exp(— A g £). For 
the gs,lp-generator the processes starting from one DO 
creating two additional DOs can be more difficult to sup- 
press if they decrease the total energy. This is possible if 
the DOs disperse and a DO at high energies decays into 
three DOs at lower energy. 

The pc-generator aims in addition at eliminating pro- 
cesses starting from two and more DOs so that there are 
even more processes which may decrease the total en- 
ergy while the number of DOs increases. Hence we are 
not surprised to see that the pc-generator induces a flow 
of the ROD which displays even more pronounced non- 
monotonic behavior. 

From these observations the conclusion to always favor 
the gs-generator suggests itself. But it is in fact a trade- 
off. The gs-generator is quicker to implement and more 
robust in its convergence, but it achieves less because it 
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decouples only the subspace without any DOs. For deriv- 
ing only an extended Heisenberg model this is completely 
sufficient and hence for this aim the gs-generator is the 
generator of choice. But if one is additionally interested 
in an explicit description of the dynamics of DOs the 
other generators are advantageous as we will illustrate 
next. 

To see how the coupling constants are influenced by the 
choice of generator a few exemplary results are shown. 
One of the most important magnetic coupling constants 
besides the nearest neighbor coupling Ji is the Heisen- 
berg interaction between next-nearest neighbors J2, i.e., 
diagonal over a plaquette. The behavior of this exchange 
coupling as function of W/U for the three generators is 



depicted in Fig. 10 



■gs 


generator 


'gs.lp 


generator 


... tl 

'pc 


generator 




FIG. 10: (Color online) Behavior of the next nearest neigh- 
bor interaction J2 for different generators as function of W/U 
determined in an upto4 calculation. The results for the three 
generators almost coincide even for large values of W/U. 

The curves for all three generators almost coincide. 
This underlines that it is completely sufficient to use the 
gs or the gs,lp-generator for the determination of this 
coupling. All processes contributing to this coupling are 
included in the CUT induced by the gs-generator. This 
can be understood from the fact, that the magnetic cou- 
pling J2 describes an interaction within the subspace of 
half-filled states, i.e., the subspace without any DO. This 
subspace is decoupled from the remainder of the Hilbert 
space by all three generators. If no truncation errors oc- 
curred, all three generators would indeed yield precisely 
the same result, cf. Ref. [T71 

In contrast to a pure spin-spin coupling the term 



H'v n - V" ^2 X! 



h.c 



} 



(27) 



acts on two DOs. It describes the hopping of an electron 
from a singly occupied site to an empty site under the 
condition that site k is occupied by a DO. It is a process 
which is not active on the subspace with zero or only 
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FIG. 11: (Color online) Density-density interaction between 
third nearest neighbors for different generators obtained in 
an upto4 calculation. The results for the pc- and the gs, lp- 
generator agree very well whereas the results of <?s-generator 
show larger deviations. 



one DO. We do not expect that the results for the dif- 
ferent generators agree. Indeed, the results for this cou- 
pling constant (Fig. [Tlj ) show rather large deviations of 
the results obtained by the gs-generator from the results 
obtained by the other two generators. This illustrates 
that the gs-generator induces a different unitary trans- 
formation than the other two generators. Note that the 
deviations do not necessarily imply that the gs-result is 
less accurate because it results from the representation 
of the Hamiltonian in a different basis. 

In view of the above arguments, it is surprising that the 
results of the gs,lp-generator and the pc-generator are so 
close to each other. From their definitions we expect that 
the pc- and the gs,lp-results agree very well for processes 
involving a single DO, but not necessarily for processes 
involving two DOs. 

In conclusion, the question which generator is opti- 
mum cannot be answered generally. It depends on the 
particular objective of the intended investigation. 



VI. APPARENT CHARGE GAP 

In the previous sections we have started to discuss the 
issue for which conditions the intended mapping from the 
Hubbard model to a generalized t-J model is possible and 
justified. Qualitatively it is obvious that U must be large 
enough. But quantitative indicators are needed. Here we 
aim at giving a quantitative estimate for the parameter 
range in which the mapping is justified. 



A. General Considerations 

The basis of the transformation from the Hubbard 
model to the t-J model is the elimination of charge fluctu- 
ations. These charge fluctuations correspond to changes 
in the number of DOs. The corresponding processes are 
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the ones that we consider to be off-diagonal. To be able 
to eliminate such processes the subspaces with differing 
numbers of DOs have to be separated in energy. In Sect. 
Ill we introduced the residual off-diagonality (ROD) as a 
measure to decide if for a given £ the off-diagonal terms 
are small enough to be neglected. 
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FIG. 12: (Color online) Behavior of the residual off- 
diagonality (ROD) for the pc-generator obtained in a double 
plaquette calculation for different values of W/U. The ROD 
measures to which extent charge fluctuations are eliminated 
as function of the flow parameter £. 

Figure [12] shows the behavior of the ROD of the pc- 
generator for the double plaquette calculation at various 
values of W/U. For small values of the flow parameter 
IU < 2.24 the ROD decreases exponentially. If the ratio 
W/U is increased to W/U = 1.12 the ROD evolves iron- 
monotonically. For £U < 64 the ROD falls below 10 
to rise again for larger values of I. For even higher values 
of W/U this increase sets in for smaller and smaller I. 

We argued in the preceding section that a non- 
monotonic behavior of the ROD as observed in Fig. [12] is 
a first clue for a possible breakdown of the mapping. Al- 
though the non-monotonicity indicates problems of the 
mapping already for W/U = 1.12 no sign of a possible 
breakdown can be seen in the coupling constants, see for 
instance Fig 
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This is due to the fact, that the dom- 
inant spin coupling constants are already converged to 
their value in the effective model for small values of £. 
Thus processes appearing only at large £ have no influ- 
ence on these values. 

To make progress in determining the range of valid- 
ity of the mapping we have to study the separation of 
energy scales between the states without DOs and the 
states with DOs. To this end we investigate the energies 
of an added electron or an added hole to the system, that 
means the density-of-states (DOS) which comprises the 
LHB and the UHB for large U, see Figs. [I] and [2j If W 
is increased the bands approach each other and eventu- 
ally touch so that there is no energy separation anymore. 
States with differing number of DOs have the same en- 
ergy so that charge fluctuations can not be eliminated. In 



a paramagnetic description by dynamic mean-held theory 
(DMFT) the Mott- insulating phase becomes unstable at 
this very pointPEMI] 

Turning the argument around we use the gap A g be- 
tween the lower and the upper Hubbard band as quan- 
titative indicator of the energy separation of subspaces 
with different number of DOs. If A g is finite there is 
good physical reason to regard the mapping of the Hub- 
bard model onto the t-J model as justified. If A g vanishes 
the mapping has to break down. 

There is one additional aspect to which we have to 
draw the reader's attention. In infinite dimensions, where 
DMFT is exact, one may suppress long-range magnetic 
order and consider the paramagnetic phase which does 
not show an spin-spin correlation between different sites 
so that the insulating paramagnetic phase behaves like 
the reference ensemble ^ for 6 — 0. In particular, the 
charge gap does not depend on the spin state. 

But in finite dimensions, even without long-range mag- 
netic order one has to expect that the charge gap will 
generically depend on the spin state. Hence there is no 
well-defined charge gap without specifying the state of 
the spin degrees of freedom. Thus we have to introduce 
the concept of the apparent charge gap A g ® which is de- 
signed to describe the energy separation of subspaces 
with different number of DOs if an electron is added 
to the disordered reference ensemble ([9]). The appar- 
ent charge gap is not rigorously defined and it cannot 
be measured precisely in experiment because it does not 
capture bandtails of the Hubbard bands which carry lit- 
tle spectral weight. But it is an estimate for the energy 
separation between states without DOs and states with 
DOs. Hence it provides an appropriate estimate of the 
range of validity of the mapping from the Hubbard model 
to the t-J model. 

We calculate A g for the effective t-J model derived 
before. In the half-filled case the gap is calculated by 
estimating the lowest possible energy of an added DO, 
for details see below. Calculations for the half-filled case 
indicate a closure of the gap for W/U « O.jP. In the 
doped case the calculation of A g is divided into two steps 
as can be understood from the DOS sketched in Fig. [2] 

To calculate the apparent gap we first determine the 
lowest possible energy Auhb of a DO in the upper Hub- 
bard band. In a second step we calculate the maximum 
energy for the destruction of a DO A l hb, that means for 
adding an electron to an empty site. Hence in the doped 
case we use 



— ^UHB — ^LHB, 

while in the undoped case we have 
A g = 2 Auhb- 



(28) 



(29) 



Note that the seemingly discontinuous definition of A g as 
function of Auhb stems from the discontinuous evolution 
of the Fermi level which jumps upon hole doping from 
the middle between the Hubbard bands to the edge of 
the lower Hubbard band. 
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B. Calculation of A 9 

The apparent charge gap is calculated for the effec- 
tive t-J model. A full diagonalization of the Hamilto- 
nian is not feasible and it would not provide what we 
need, namely the charge gap above the disordered refer- 
ence ensemble ([9]). Hence we apply a Lanczos approach 
in terms of operators. The Lanczos approach is appropri- 
ate because we only aim at extremum eigenvalues. Since 
we have to deal with operators acting on the reference 
ensemble^, which is a mixed state, the Liouville formu- 
lation of this method has to be usecPsHUl. The evolution 
of an operator A is given by the Liouville superoperator 



LA = 



H r ff, A 



(30) 



The effect of this superoperator applied to the creation 
operator of a DO consists of moving the DO and chang- 
ing the spin background. With this operator a basis of 
operators {vo, ...v n } describing the DO with momentum 
k and the effect on its surrounding spins is built recur- 
sively. In the hrst part of the calculation the minimal 
energy of a DO with momentum k is calculated. The 
calculation starts with the vector 



1 



N H 



E 



e ikp fi r t 



(31) 



where N denotes the number of lattice sites and the vec- 
tor r determines the actual position of the DO. The ac- 
tion of this operator is to put an f electron on a site 
which is already occupied by a J, electron so that a dou- 
bly occupied site is created. From this starting vector 
the basis is built recursively by 



(32) 



according to the rules of the Lanczos tridiagonalization. 
The scalar product of the Liouville formulation^ is de- 
fined as 



(A\B) = Tr (&B po) . 



(33) 



with the statistical operator po of the reference ensemble 
([9]). The prefactors <Xj are given through the projection 
onto Vj 



and the bi are given by 



{Vi\£Vj) 
(Vi\Vi) 



(Vi\Vj) 
(Vi-tlVi-t) 



(34) 



(35) 



In this operator basis the Liouville superoperator takes 
tridiagonal form with the coefficients <ij on the diagonal 
the coefficients bi on the secondary diagonals. 



If infinitely many iterations were performed the dis- 
persion of a DO relative to the disordered spin back- 
ground would be given as the lowest energy in the sub- 
space spanned by the calculated operator basis. In real 
calculations only a few iterations are feasible due to the 
humongous number of terms in the effective Hamiltonian. 
Starting with the vector w in ( 31 ) consisting of one single 



operator, the commutation leads to increasingly compli- 
cated terms whose appropriate superposition describes 
Vj. The effort grows exponentially with the number of 
iterations. Thus we restrict ourselves to a finite basis 
{vq . . . v n }. The lowest energy calculated in the subspace 
spanned by {v\ . . .v n } yields an upper bound Auhb to 
the real dispersion of the DO. Note that in the following 
we denote by Auhb this upper bound in order to keep 
the notation simple. 

In the half-filled case the apparent charge gap results 
from Eq. ( 29 1 . For finite doping the particle- hole sym- 



metry is lost, the Fermi energy jumps to the LHB, and 
the value Alhb has to be determined. To this end, we 
start from the modified vector 



1 



u = 



N ^ 



( 



1 



(36) 



This operator destroys a DO by placing a single electron 
on an empty site. The value for Alhb we obtain from 
the calculation in a finite subspace {{to . . . u n } is a lower 
bound to the true maximum energy. Finally the apparent 
charge gap is given by ( 28 1 in dependence on the doping 



level S. As argued before the mapping to the t-J model 
is justified as along as A g > 0. 

Extremum values of Ek occur at the high symmetry 
points of the Brillouin zone. Thus we avoid costly calcu- 
lation of the whole dispersion and focus on the momenta 
k = (0, 0) and k = (n, ir) where the lattice constant is set 
to unity. The calculations rely on the the nearest neigh- 
bor effective model. Previous calculations in the half- 
filled case^ showed that there is no significant change in 
the results obtained for different truncation schemes be- 
cause the main uncertainty results from the limited num- 
ber of iterations in the Lanczos tridiagonalization. The 
truncation of the effective model used plays only a minor 
role. Since only a few iterations were feasible we ad- 
ditionally perform an extrapolation. The results for the 
gaps are extrapolated in 1 jn with n denoting the number 
of iterations. By extrapolating to n = 00 we obtain an 
estimate for A g ; for more details we refer to Ref. El Fig- 
ure [13] displays the extrapolation for the half-filled case. 
For the doped case, Alhb an d Auhb are extrapolated 
separately in 1/n. 



C. Results for A 9 

The apparent charge gap is computed for the effec- 
tive t-J model derived by a CUT with NN truncation 
using the pc-generator or the gs,lp-generator. The gs- 



13 




FIG. 13: (Color online) Linear extrapolation in ^ of the ap- 
parent charge gap at half-filling for W/U = 0.5 obtained by 
n iterations. 



generator is not used in this context because the result- 
ing effective model mixes a single DO with the subspaces 
of two and more DOs. The gap is calculated for various 
doping levels as function of W . Thus the value W/U 
up to which the mapping is justified is estimated from 
A g (W/U) = 0. The minimum Atjhb of the dispersion 
of a DO is found for a vanishing momentum. The maxi- 
mum energy Alhb for the destruction of a hole is found 
for a momentum k = (n, it). 




FIG. 14: (Color online) The extrapolated apparent charge 
gap as function of W/U for various 8. 



The results for the extrapolated A g are displayed in 
Fig. [14] for various values of the doping 6. For vanish- 
ing bandwidth W = the apparent gap is given by the 
Hubbard repulsion U. Thus the curves of A app /J7 start 
at unity. Then the gap decreases almost linearly until 
A app = is reached. Negative values of the gap indicate 
the breakdown of the mapping. The linear decrease of 
the gap has also been observed for the half-filled case^l. 



The decrease leads to a closure of the charge gap for 
W/U — 1. For the Bethe lattice with z->ooa closure of 
the gap was found for W/U = 0.8tP3 which agrees well 
with our estimate in view of the different lattices and 
techniques. Other numerical evaluations of DMFT for 
the Bethe lattice yield a closure of the insulating gap at 
W/U « 0.8423 or at W/U sa 0.8SP. 
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FIG. 15: (Color online) Limiting values for W/U up to which 
a mapping of the Hubbard model to the effective t-J model 
is justified. It is derived from the zeros A g (W/U) = as 
function of doping. 

Our results indicate that the apparent charge gap A g 
for the square lattice closes at W/U ps 0.98 for S = 0. 
Upon doping A g vanishes even faster upon increasing 
bandwidth so that the range of applicability of the map- 
ping 'Hubbard model — > t-J model' is reduced. From 
the values of W/U where A ff beomes zero we estimate 



this range of applicability. The result is shown in Fig. 15 
which represents one of the central results of this work. 
Our approach provides the first systematic derivation of 
this diagram of applicability as function of doping. 

The range of applicability decreases from W/U « 0.98 
for 5 = to W/U fa 0.67 for S = 0.25. Then the range 
of applicability increases again slightly to W/U = 0.73 
for S — 0.4. The plateau in W/U {8) and the moder- 
ate increase are rather unexpected, cf. RefW We do 
not have an obvious explanation for it. In contrast, the 
decrease of the range of applicability for S < 0.3 meets 
the qualitative expectation since a doped system has a 
higher mobility of charges so that the energy separation 
of sectors of differing number of DOs becomes smeared 
out. 

The relative constant limiting value for W/U below 
which the use of a generalized t-J model is justified pro- 
vides interesting information on the applicability of t-J 
models for doped systems. The use of t-J models is very 
widespread in theoretical studies for the high-T c super- 
conductors based on cuprates. Commonly used parame- 
ters are W/U ps 0.7 and S < 0.«P^. Our results indicate 
that the use of t-J models is indeed justified. But caution 
is required in the doping range 0.18 j§ S ;$ 0.25 where 
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W/U « 0.7 is at about the limit of applicability. Thus 
our results shed light on the important question of the 
applicability of a commonly used model. It is remarkable 
that the issue of how to justify this model for significant 
levels of dopings has attracted so little attention so far. 



VII. RESULTS FOR THE RELEVANT 
COUPLING CONSTANTS 



In the preceding section we comprehensively discussed 
the applicability of the derivation of a generalized t-J 
model. The result of this discussion is summarized in 
the estimated range of applicability shown in Fig. 15 In 



the present section, we provide the coupling constants 
which ensue from the CUT of the Hubbard model to the 
t-J model. Results are given in the range W/U < 1.0 
because the mapping definitely breaks down beyond. 

All results shown are derived from upto4 calcula- 
tions using the pc-generator. Additionally we performed 
random double plaquette calculations with the gs,lp- 
generator to check if there are changes in the coupling 
constants when higher truncation schemes are applied. 
No significant differences are found. Thus the upto4 trun- 
cation appears to be sufficient to determine the coupling 
constants. The results for the half-filled case (<5 = 0) 
shown in the following figures agree perfectly with the 
results obtained by Reischl et alP. 

Although the mapping generates a large number of 
terms, only few of them are really relevant in the final 
effective model. Most others only have very small pref- 
actors. 



Chemical Potential 



First, we consider the chemical potential as defined in 
(24). In leading quadratic order of jj it is proportional to 



(2) 

5. Therefore /i will be shown in units of 6J{ as defined in 



( 12 1. The ratio fi/5 shows almost no dependence on S, see 
Fig. 16 As function of W/U the chemical potential stays 
rather constant and even for W/U = 1 the deviations are 
small. The dependence on W/U is greater for smaller 
values of <5 concentration. 




FIG. 16: (Color online) Chemical potential relative to <5j{ 2 'for 
various doping concentrations 8. 



The left panel of Fig 
the ratio W/U 



17 shows the dependence of J\ 
Starting from 



on tne ratio W/U. Starting irom J{ ' for W/U = 
the coupling constant takes slightly smaller values for 
larger W/U. Additionally, the doping dependence of J\ 
for various values of W/U is shown relative to its value 
for the undoped system in the right panel of Fig. |17| J\ 
increases with S. The doping has a greater influence for 
larger values of W/U, but the effect remains rather small. 
Even for W/U — 0.8 the doping causes a change in Ji of 
only about 3%. 




FIG. 17: (Color online) Dependence of J\ on W/U for various 
values of <5 (left panel). The dependence on S for two ratios 
W/U is depicted in the right panel. The undoped value for 



W/U = 0.05 is found to be Ji(0) « 1.5621 
W/U = 0.8 Ji(0) takes the value 0.037917. 



W 1 U. For 



B. Spin Terms 

The dominant terms of the effective model are the 
Heisenberg-type spin interactions. 

HlleiBenb.\i-j\ = J\i-j\SjSj ■ (37) 

The largest contribution of this type is the Heisenberg 
exchange J\ between nearest neighbors. All results are 

(2) 

shown relative to the leading perturbative result J{ = 

4ti 
u • 



The Heisenberg exchange between next-nearest (diag- 
onal) neighbors J2 as well as the exchange J3 between 
neighbors at a linear distance of two lattice spacings are 
much smaller than J\ . Both terms appear in fourth order 
of jj. Even for W/U = 1 J2 is smaller than 0.03Ji, see 
Fig. |18| Surprisingly, J2 shows a slightly more significant 
(relative) dependence on the dopant concentration than 
J\, see Fig. [18] Yet, in view of the small absolute values 
of J2, this doping dependence can be neglected. 

As can be seen in the right panel of Fig. [l9]the coupling 
J3 shows a counter-intuitive behavior. First it decreases 
upon doping but increases again beyond 5 « 0.3. 
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FIG. 20: (Color online) Positions of spins interacting via the 
ring exchange Hq. 



in an effective model. 



The ring exchange shown in Fig. 21 displays nearly no 
dependence on the doping S. Even for doping as large as 
6 = 0.8 the change in the coefficient is less than 1.12 per- 
cent. Thus while ring exchange is an important process 
its doping dependence can safely be omitted. 



FIG. 18: Effective J 2 for various doping concentrations as 
function of W/U. 




FIG. 19: (Color online) Dependence of J3 on W/U for various 
values of 8 (left panel). The dependence on 5 for two ratios 
W/U is depicted in the right panel. The undoped values are 
J 3 (0) « 6.1039 • 10" 9 U for W/U = 0.05 and J 3 (0) w 3.9802 • 
10" 4 U for W/U = 0.8. 



Besides the 2-spin terms in ( 37 1 the two dimensional 

The 



square lattice also allows for 4-spin interactions 
leading contribution is given by 



<i,j,k,l> 



(38) 



which we sloppily call ring exchange although the com- 
plete ring exchange comprises also nearest-neighbor and 
diagonal 2-spin couplings^. We do so since these 2-spin 
terms are accounted for by J\, J2, and J3 in ([37]). The 



term (|38|) describe s th e interaction of the four spins on a 
and it occurs first in order (jj) 20. 
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plaquette, see Fig. 

Its importance is discussed at length in the literature, 
see for instance. Refs. 9 34-38 and references therein. 
The magnetic excitations in planar cuprates may not be 
understood without considering ring exchange^HlU. 

Compared to other quartic exchange couplings such as 
J2 or J3 the ring exchange is much more important, see 
its values in Fig. [21] The ring exchange takes values of 
up to 20% of J\. Hence this term must not be neglected 




FIG. 21: (Color online) Effective ring exchange term Jo for 
various values of S (left panel) and its doping dependence of 
Ju (right panel) relative to the undoped values Jp ~ 2.4390 • 
10 ( - 7) U for W/U = 0.05 and J n (0) « 0.0125 U for W/U = 
0.8. 

The second 4-spin term is the cross exchange 

H x = J x {^Sk) (4$) • ( 39 ) 

<i,j,k,l> 

In this term the spins are located on the same sites as for 
the ring exchange, but the inner products are taken of 
the diagonal spins. The corresponding coupling constant 
is shown in Fig. [22] It takes smaller values than Jq and 
hardly shows any doping dependence. 

C. Interaction of Double Occupancies 

The effective generalized t-J model also contains in- 
teractions between DOs. First, we consider the Hubbard 
repulsion U which determines the energy costs for the 
creation of a single DO. So strictly speaking it does not 
represent a true interaction. Since the deviations of the 
doped values of U from the ones in the half-filled case 
are small we directly show the doped values relative to 
the half-filled ones in Fig. [23] This coupling constant 
shows nearly no dependence on the doping S. Hence the 
influence of doping on U may be neglected. 

The following interaction terms are active only in the 
presence of at least two DOs. Thus these terms have 
to be seen as genuine 2-DO interactions. Among them 
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FIG. 22: (Color online) Effective J x as denned in Eq.|39]for 
various doping concentrations 5. 
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FIG. 23: (Color online) Hubbard repulsion for the doped case 
relative to the value in the half-filled case as function of 5 for 
various values of W/U. The doped values are given relative 
to the undoped values of [7(0) « 1.0003(7 for W/U = 0.05 
and [7(0) w 1.078717 for W/[/ = 0.8. 



density-density interactions of various distances appear. 
The density-density interaction between nearest neigh- 
bors reads 



V ^2 TH.snjj, 

<i,j> 



(40) 



where hs denotes the operator counting the number of 
electrons on a site compared to the average filling, cf. 
Tab. HJ At half- filling this term only contributes if site 
i and site j are either empty or truly doubly occupied. 
Figure [24] only shows an increase in the coupling constant 
of about 1% under the influence of doping for W/U = 0.8. 

A second type of interaction is correlated hopping. 
The most important term of this kind is the hopping of 
two electrons to a nearest neighbor site which is initially 



FIG. 24: (Color online) Effective density-density interaction 
as function of doping for various values of W/U. Results are 
given relative to the values at zero doping V(0) w —3.9056 • 
10~ 5 U for W/U = 0.05 and V(0) w -9.2525 • 10~ 3 U for 
W/U = 0.8. 



empty, see Eq. 22 Since the empty state also corresponds 
to a DO, the effect of the term is to exchange the two 
DOs. The results for various doping levels are depicted 
in Fig. [25] 




FIG. 25: (Color online) Effective pair interaction V v as defined 
in Eq. 
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for various 5 as a function of W/U . 



Besides the nearest neighbor pair interaction V p there 
are also pair interaction terms between three spins. One 
of these terms is the interaction of three spins on a pla- 
quette which reads 



<? <i,j,k> 



(41) 



The sites i and j are supposed to be diagonal neighbors 
with a common adjacent site k. One possible process con- 
sists of the hopping of an electron from a singly occupied 
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site j to an empty site k. Simultaneously, an electron 
from the doubly occupied site i hops to site k forming 
a DO on this site. The corresponding effective coupling 
constant Vp is depicted in Fig. 26 as function of doping. 



Since this correlated hopping imposes an additional 
constraint on the state of site k it is half as large as the 
nearest neighbor term V p . Even for large values of W/U 
the coupling constant is increased only by 8% for large 
doping. 





... w/u 


= 0.05 




-- w/u 


= 0.8 



X 1.04 
> 1.02 



FIG. 26: (Color online) Effective pair interaction Vp on W/U 
for various 5. In the undoped case Vp takes the value Fp'(0) « 
3.9035 ■ 1CT 5 U for W/U = 0.05 and F p '(0) w 8.5167 ■ 10~ 3 U 
for W/U = 0.8. 

The last class of terms considered are correlated hop- 
ping terms such as 



\ a A*\fl>j,p fi k + h - c -} • ( 42 ) 

One of the processes described by Hy' is the hopping of 
an electron from a singly occupied site j to an empty site 
i under the condition that site k is occupied by a DO, 
see Fig. [27j Sites i and j are diagonal neighbors on a 
plaquette and k joint adjacent neighbor. 



FIG. 27: (Color online) Example for processes comprised in 

Due to the constraint that site k has to be occupied by 
a DO and site i has to be empty in the beginning, these 
processes rely on the presence of two DOs which justifies 
to view them as true interactions. The number of DOs 



is not changed by this process. Processes such as Hy 
appear in second order of jj. 

The corresponding coupling constant is shown in the 
left panel of Fig. 28 as function of W/U . In the right 
panel of Fig. 



2.s 



the value for the coupling constant V„ 
in the doped case is shown relative to its value in the 
half-filled case. For large values of W/U, the coupling V„ 
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FIG. 28: (Color online) Effective correlated hopping process 
V^ over diagonal neighbors for various 8 (left panel) and as 
function of 5 for the values W/U = 0.05 and W/U = 0.8 
(right panel). The undoped values are V„ ~ -1.9518 ■ 10" 5 U 
(W/U = 0.05) and V/> « -4.2584 ■ 10" 3 U (W/U = 0.8). 

shows a noticeable dependence on S. Note that besides 



the correlated hopping defined in (|42|) and illustrated in 
Fig. 
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there is correlated hopping between three sites 
located on three sites in a row. The corresponding cou- 
pling constant V/[ shows the same behavior as so that 
we do not show it here for brevity. 

For too large values of W/U, i.e., W/U ^> 1, the curves 
for the coupling constants are not smooth anymore (not 
shown here^P. This observation is explained by the break- 
down of the mapping as it is indicated by the behavior 
of the apparent charge gap (see Sect. VI I. 



D. Hopping Terms 

The first term to be considered is To which was intro- 
The initial Tq represents hopping 



duced before in (7a 



processes by one lattice spacing without a change in the 
number of DOs. The corresponding coupling constant to 
is shown in the left panel of Fig.[29]relative to its unrenor- 
malized value. The deviation of the coupling to from its 
bare value is proportional to (^) . It is increased upon 
increased doping. To examine the doping dependence the 
renormalized value of to i n the doped case is compared 
to the one in the half-filled case in the right panel of Fig. 



29 Under the influence of doping the hopping param- 
eter is increased linearly. But even for W/U = 0.8 the 
parameter is changed only by a few percent. 

Hopping also occurs between diagonal sites on a pla- 
quette, for instance in 



rpl 



t! 



E E 

a <<i,j>> 



[(' 



c . - C - 



,(1 



(43) 
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is 



where the sum runs over two diagonal neighbors i and 
j which have a common nearest neighbor k. The size of 
the corresponding hopping parameter t' spin , see Fig, 
comparable to the spin independent parameter t'. 
shows that the induced NNN spin dependent hopping 
processes are as important as the spin independent ones. 
This was first observed by Reischl et alP at half-filling. 



FIG. 29: (Color online) Effective NN hopping without chang- 
ing the number of DOs compared to its initial, bare value 
(left panel). Its doping dependence relative to the values 
at half-filling t (0) ~ 6.2490 • 1(T 3 U for W/U = 0.05 and 
t (0) » 0.0956 U for W/U = 0.8 is depicted in the right panel. 



Here the double bracket under the sum indicates next- 
nearest neighbors (NNN) on the square lattice. The 
same process also appears between third nearest neigh- 
bors with a distance of two lattice sites. The correspond- 
ing coupling constant is denoted by t" . Since t" and t' 
show very similar behavior we only show the results for 
t' in Fig. [30 
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FIG. 30: (Color online) Effective NNN Hopping t' between 
diagonal sites as function of W/U for various values of 8 
(left panel) and its doping dependence for W/U = 0.05 and 
W/U — 0.8 (right panel). In the undoped case t' takes 
the values i'(0) « -3.9053 ■ 10" 5 U (W/U = 0.05) and 
t'(0) w -9.4785 • 10" 3 U for W/U = 0.8. 



The hopping t' decreases linearly for increasing ratio 
W/U with slopes depending on the doping, see left panel 
of Fig. [30] Relative to its values at half-filling the de- 
crease as function of doping hardly depends on W/U, see 
right panel of Fig. 30 It is remarkable that the con- 
stant is decreased to almost for 6 — » 1. This is actually 



the only significant dependence on doping that we found. 
But one has to keep in mind that the absolute values of 
t' are small. Note that the sign of to is positive whereas 
t' and t" are negative. 

An interesting coupling generated in second order of 
if is the spin dependent hopping described by 



f ■ = t' ■ V V 

spin spin / j / j 



a/3 <i,fe,j> 



[[ 



1 - n 



i, a J ^l&^c^il-n^) 
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FIG. 31: (Color online) Effective spin dependent hopping t' spin 
for various values of 8 (left panel) and its dependence on 8 for 
two values of W/U, namely W/U = 0.05 {t' spin (0) « 3.9042 ■ 



10~ 5 U) and W/U = 0.8 (t' sf 
panel) . 



, n (0) « 8.8496 • W' s U) (right 



S, 



(44) 



Upon doping the spin dependent hopping term is in- 
creased. But even for larger values of W/U the value of 
t' spin is increased only by a few percent. The analogous 
process also appears between third nearest neighbors. 
The corresponding coupling constant behaves similar to 
f^pjj, so that we do not display it here. Both processes 
concern three sites. Thus it does not matter significantly 
whether the sites are aligned linearly or in a right angle 
on a plaquette. 

Compared to the spin independent hopping t' the spin 
dependent hoppings t' spin and £" pin are not negligible. But 
all of them are fairly small compared to the bare NN hop- 
ing to- Thus one can either stick to a pure t-J model or 
include more extended hopping terms. But if one opts for 
including hopping over two lattice spacings one should in- 
corporate spin independent hopping t' as well as the spin 
dependent hoppings t' spin and t" pin . The doping depen- 
dence of the spin dependent hopping elements may be 
neglected. In contrast, the hopping element t' shows a 
rather strong dependence on the doping concentration 5. 



VIII. SUMMARY 

We presented a systematically controlled mapping of 
a fermionic Hubbard model to a generalized t-J model. 
The conceptual foundation of this mapping was analyzed 
carefully. In particular, we developed a scheme for this 
mapping which covers also the interesting case of sub- 
stantial doping. Remarkably, this issue had so far at- 
tracted only little attention. 

In the derivation of the generalized t-J model we elim- 
inate the charge fluctuations by self-similar continuous 
unitary transformations. Processes that change the num- 
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ber of double occupancies are rotated away. Thereby, we 
obtain the effective coupling constants as function of the 
doping 5 and of the ratio W/U where W is the bandwidth 
and U the local interaction. Note that the generalized t- 
J model comprises the magnetic degrees of freedom as 
well as the kinetics and the interactions of double occu- 
pancies. 

We extended the concept of the apparent charge gap 
AjSl from half-filling to the doped system. This gap is 
not the true physical gap but it measures the energy 
separation of subspaces with differing number of dou- 
ble occupancies irrespective of the spin state. We argue 
that as long as A g is finite the mapping to a t-J model 
is justified. A vanishing A g indicates the breakdown 
of this mapping. By estimating the parameter where 
A g (W/U, 5) — holds we derived a diagram of applica- 



bility of the t-J model shown in Fig. 15 As expected the 
applicability is reduced upon doping 5. But it levels at 
intermediate values of doping so that the commonly as- 
sumed parameters for the description of high-T c cuprates 
lie within the range of applicability. To our knowledge, 
no such result was derived before. 

Furthermore we find that the coupling constants of the 
effective model show hardly any doping dependence. The 
only coupling which exhibits a significant dependence on 
6 is the hopping parameter t' describing hopping between 
diagonal neighbors. Relative to its value at half-filling t' 
exhibits a strong doping dependence. But the absolute 
value of this hopping element remains small. Thus within 
a wide range of doping the t-J model with constant cou- 
pling constants is appropriate. Besides the usually con- 
sidered terms, the 4-spin ring exchange on each plaquette 
should be included. 

Technically we used recently developed types of 
infinitesimal generators for the continuous unitary 
transformational. They only decouple certain subspaces 
of the Hilbert space which simplifies and accelerates the 
calculations. So far, the pc-generator was used which 
leads to a particle number conserving effec tive model; the 
particles are the double occupancie^^l. We extended 
the gs-generator introduced previously for the ground 
state of matrice^H and of many-body systems^ to mixed 
reference ensembles. The gs-generator is particularly 
suited to obtain the purely magnetic Heisenberg model 
since it efficiently decouples the subspace of the refer- 
ence ensemble from the remainder of the Hilbert space. 
If, however, the dynamics of the double occupancies mat- 
ters as well, the gs,lp-generator turned out to be a good 
compromise between efficiency and sufficient decoupling. 



This generator decouples the reference ensemble and the 
states with one double occupancy from the rest of the 
Hilbert space. We found that the couplings derived from 
a faster gs,lp-calculation agree very well with the results 
from a slower pc-calculation. We expect that these gen- 
erator or modifications of them will continue to play an 
important role in the systematic derivation of effective 
models. 

The present analysis for the square lattice can certainly 
be extended to other types of lattices such as the trian- 
gular lattice which has already been analysed by per- 
turbative CUT^H, the honeycomb lattice^, or more so- 
phisticated lattices such as the kagome lattice and so on. 
In this way, the effects of subleading magnetic exchange 
processes such as ring exchange can be analysed quanti- 
tatively. 

Another route to extend the present calculation is to 
also transform the observables, for instance the standard 
fermionic creation operator. At half-filling, one will then 
be able to compute the spectral weight in the upper Hub- 
bard band, that means in the subspace with one double 
occupancy. But there should be also weight in the sub- 
spaces with three and more double occupancies. To our 
knowledge, no estimate whatsoever exists for the weight 
in such trans-Hubbard bands. 

More generally, the systematic derivation of effective 
models in other contexts can be tackled by adapting 
the ideas of the present work. For instance, the reliable 
downfolding of interacting fermionic models with many 
bands to models with a minimum number of bands and 
Hubbard-type of interactions is a long standing field of 
researdPSHSl w hich continues to attract much attention, 
see for instance Refs. I47H501 We think that continuous 
unitary transformations provide an promising approch to 
make systematic and controlled progress in this field. 

Hence we expect that the systematic derivation of ef- 
fective models by means of continuous unitary trans- 
formations will continue to evolve into a field with 
widespread applications. 
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